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Abstract. With the aim of understanding the emergence of collective motion from lo- 
cal interactions of organisms in a "noisy" environment, we study biologically inspired, 
0^ ■ inherently non-equilibrium models consisting of self-propelled particles. In these models 

particles interact with their neighbors by turning towards the local average direction of 

' (-H ' ] motion. In the limit of vanishing velocities this behavior results in a dynamics analogous 

O , . to some Monte Carlo realization of equilibrium ferromagnets. However, numerical sim- 

I ' ulations indicate the existence of new types of phase transitions which are not present in 

^ , the corresponding ferromagnets. In particular, here we demonstrate both numerically 

(^ ■ and analytically that even in certain one dimensional self-propelled particle systems 

J ' an ordered phase exists for finite noise levels. 

o ' 

^' 1 Introduction 

^ ; The collective motion of organisms (flocking of birds, for example), is a fascinat- 

ing phenomenon capturing our eyes when we observe our natural environment. 
In addition to the aesthetic aspects, studies on collective motion can have inter- 
^ , esting applications as well: a better understanding of the swimming patterns of 

C^ ' large schools of fish can be useful in the context of large scale fishing strategies, 

Cn , or modeling the motion of a crowd of people can help urban designers. Here we 

address the question whether there are some global, perhaps universal features 
of this type of behavior when many organisms are involved and such parame- 
^^ ' ters as the level of perturbations or the mean distance between the organisms is 

Q^ i changed. 

^^ ' Our interest is also motivated by the recent developments in areas related 

O i to statistical physics. During the last 15 years or so there has been an increas- 

C/3 ' ing interest in the studies of far-from-equilibrium systems typical in our natural 

^^1 and social environment. Concepts originated from the physics of phase transi- 

tions in equilibrium systems [1] such as collective behavior, scale invariance and 
renormalization have been shown to be useful in the understanding of various 
non-equilibrium systems as well. Simple algorithmic models have been helpful 
in the extraction of the basic properties of various far-from-equilibrium phe- 
^ • nomena, like diffusion limited growth [2], self-organized criticality [3] or surface 

5^ I roughening [4]. Motion and related transport phenomena represent a further 

characteristic aspect of non-equilibrium processes, including traffic models [5], 
thermal ratchets [6] or driven granular materials [7] . 

Self-propulsion is an essential feature of most living systems. In addition, 
the motion of the organisms is usually controlled by interactions with other 
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organisms in their neighborhood and randomness also plays an important role. In 
Ref. [8] a simple model of self propelled particles (SPP) was introduced capturing 
these features with a view toward modeling the collective motion [9, 10] of large 
groups of organisms such as schools of fish, herds of quadrupeds, flocks of birds, 
or groups of migrating bacteria [11, 12, 13, 14, 15], correlated motion of ants [16] 
or pedestrians [17]. Our original SPP model represents a statistical physics-like 
approach to collective biological motion complementing models which take into 
account much more details of the actual behavior of the organism [9, 18], but, 
as a consequence, treat only a moderate number of organisms and concentrate 
less on the large scale behavior. 

In spite of the analogies with ferromagnetic models, the general behavior of 
SSP systems are quite different from those observed in equilibrium models. In 
particular, in the case of equilibrium ferromagnets possessing continuous rota- 
tional symmetry the ordered phase is destroyed at finite temperatures in two 
dimensions [19]. However, in the 2d version of the SSP model phase transitions 
can exist at finite noise levels (temperatures) as it was demonstrated by simula- 
tions [8, 20] and by a theory of flocking developed by Toner and Tu [21] based on 
a continuum equation proposed by them. Further studies showed that modeling 
collective motion leads to similar interesting specific results in all of the relevant 
dimensions (from 1 to 3). Therefore, after introducing the basic version of the 
model (in 2d) we discuss the results for each dimension separately and then focus 
on the Id case which is better accessible for theoretical analysis. 

2 A generic model of two dimensional SPP system 

The model consists of particles moving on a plane with periodic boundary condi- 
tion. The particles are characterized by their (off-lattice) location x^ and velocity 
Vi pointing in the direction ??i. To account for the self-propelled nature of the 
particles the magnitude of the velocity is fixed to wq . A simple local interaction 
is defined in the model: at each time step a given particle assumes the average 
direction of motion of the particles in its local neighborhood S{i) with some 
uncertainty, as described by 

Mt + ^t) = {m)s(^)+C, (1) 

where the noise ^ is a random variable with a uniform distribution in the interval 
[—rj/2,ri/2]. The locations of the particles are updated as 

x,(t + Z\t)=x,(i)+v,(t)Z\t. (2) 

The model defined by Eqs. (1) and (2) is a transport related, non-equilibrium 
analog of the ferromagnetic models [22]. The analogy is as follows: the Hamil- 
tonian tending to align the spins in the same direction in the case of equilibrium 
ferromagnets is replaced by the rule of aligning the direction of motion of par- 
ticles, and the amplitude of the random perturbations can be considered pro- 
portional to the temperature for rj <S^ 1. From a hydrodynamical point of view, 
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in SPP systems the momentuni of the particles is not conserved. Thus, the flow 
field emerging in these models can considerably differ from the usual behavior 
of fluids. 



3 Collective motion 

The model defined through Eqs. (f) and (2) was studied by performing large- 
scale Monte-Carlo simulations in Ref. [20]. Due to the simplicity of the model, 
only two control parameter should be distinguished: the (average) density of 
particles g and the amplitude of the noise 77. 

For the statistical characterization of the configurations, a well-suited or- 
der parameter is the magnitude of the average momentum of the system: <f> = 

^ ■ Vj /TV. This measure of the net flow is non-zero in the ordered phase, and 

vanishes (for an infinite system) in the disordered phase. 

Since the simulations were started from a disordered configuration, (/)(t = 
0) « 0. After some relaxation time a steady state emerges indicated, e.g., by the 
convergence of the cumulative average (1/r) J„ (t){t)dt. The stationary values of 
(j) are plotted in Fig. la vs ry for g = 2 and various system sizes L. For weak noise 
the model displays long-range ordered motion (up to the actual system size L, 
see Fig. 2), that disappears in a continuous manner by increasing 77. 

As L ^ 00, the numerical results show the presence of a kinetic phase tran- 
sition described by 

^{,)^H'^y for, <,.(,) (3) 

I for ?7 > ric{g) 

where ?7c(£') is the critical noise amplitude that separates the ordered and disor- 
dered phases and (i ~ 0.42 ± 0.03, found to be different from the the mean-field 
value 1/2 (Fig lb). 

In an analogy with equilibrium phase transitions, singular behavior of the or- 
der parameter as a function of the density and critical scaling of the fluctuations 
of (j) was also observed. These results can be summarized as 

<i){ii,g) = 4>{vl'nc{Q)), (4) 

where (i>{x) ~ (1 — x)^ for x < 1, and ^(x) « for x > 1. The critical line rjc{g) 
in the rj — g parameter space was found to follow 

vM - g", (5) 

with K = 0.45 ± 0.05. Eq.(4) also imphes that the exponent /3', defined as (ji ~ 
{q ~ Qc)^ for Q > Qc [8], must be equal to /3. The standard deviation (a) of the 
order parameter behaved as 

cr(?7) ~ |?7-r7cr'^, (6) 
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Fig. 1. (a) The average momentum of the 2D SPP model in the steady state vs the 
noise amplitude rj for g — 2 and four different system sizes [{<>) N = 800, L = 20; (+) 
N = 3200, L = 40; (□) iV = 20000, L = 100 and (x) A^ = 10^ L = 223]. (b) The 
order present at small rj disappears in a continuous manner reminiscent of second order 
phase transitions: (f> ~ [{ric{L) — ??)/»?c(i)] = (Arj/ric) , with (3 = 0.42, different from 
the mean- field value 1/2 (solid line). 
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Fig. 2. The velocities of the SPPs are displayed for various values of density and noise. 
The actual velocity of a particle is indicated by a small arrow, while its trajectory for 
the last 20 time step is shown by a short continuous curve. For comparison, the range 
of the interaction is displayed as a bar. (a) At high densities and small noise {N — 300, 
L = 5 and rj = 0.1) the motion becomes ordered, (b) For small densities and noise 
{N = 300, L = 25 and 77 = 0.1) the particles tend to form groups moving coherently 
in random directions. 
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with an exponent 7 close to 2, which value is, again, different from the mean-field 
result. 

These findings indicate that the SPP systems can be quite well characterized 
using the framework of classical critical phenomena, but also show surprising 
features when compared to the analogous equilibrium systems. The velocity vq 
provides a control parameter which switches between the SPP behavior {vq > 0) 
and an XY fcrromagnct (uo = 0). Indeed, for wo = Kosterlitz-Thouless vortices 
[23] could be observed in the system, which turned out to be unstable for any 
nonzero vq investigated in [20]. 

4 Phase diagram for a 3d SPP system 

In two dimensions an effective long range interaction can build up because the 
migrating particles have a condiderably higher chance to get close to each other 
and interact than in three dimensions (where, as is well known, random trajec- 
tories do not overlap). The less interaction acts against ordering. On the other 
hand, in three dimensions even regular ferromagnets order. Thus, it is interesting 
to see how these two competing features change the behavior of 3d SPP systems. 
The convenient generalization of Eq. (1) for the 3d case can be the following [24]: 

v,(t + At) = vo N( N((v(t))s(,)) + C), (7) 

where N(u) = u/|u| and the noise ^ is uniformly distributed in a sphere of 
radius 77. 

Generally, the behavior of the system was found [24] to be similar to that 
of described in the previous section. The long-range ordered phase was present 
for any g, but for a fixed value of rj, (j) vanished with decreasing g. To compare 
this behavior to the corresponding diluted ferromagnet, 4>{r], g) was determined 
for Wo = 0, when the model reduces to an equilibrium system of randomly dis- 
tributed "spins" with a ferromagnetic- like interaction. Again, a major difference 
was found between the SPP and the equihbrium models (Fig. 3): in the static 
case the system does not order for densities below a critical value close to 1 which 
corresponds to the percolation threshold of randomly distributed spheres in 3d. 

5 Ordered motion for finite noise in Id 

In order to study the Id SPP model a few changes in the updating rules had to 
be introduced. Since in Id the particles cannot get around each other, some of 
the interesting features of the dynamics arc lost (and the original version would 
become a trivial cellular automaton type model). However, if the algorithm is 
modified to take into account the specific crowding effects typical for Id (the 
particles can slow down before changing direction and dense regions may be 
built up of momentarily oppositely moving particles) the model becomes more 
realistic. 
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Fig. 3. Phase diagram of tlie 3d SPP and tiic corresponding ferromagnetic system. 
The diamonds show our estimates for the critical noise for a given density for the SPP 
model and the crosses show the same for the static case. The SPP system becomes 
ordered in the whole region below the curved line connecting the diamonds, while in 
the static case the ordered region extends only down to the percolation threshold p ~ 1. 



Thus, in [25] N ofF-lattice particles along a line of length L has been consid- 
ered. The particles arc characterized by their coordinate Xi and dimcnsionless 
velocity u^ updated as 



Xi{t + At) = Xi{t) + VQUi{t)At, 

U,{t + At)=G({u{t))si,))+^^. 



(8) 
(9) 



The local average velocity {u)gu\ for the ith particle is calculated over the par- 
ticles located in the interval [xi — A,Xi + A], where we fix Z\ = 1. The function 
G incorporates both the propulsion and friction forces which set the velocity in 
average to a prescribed value vq: G{u) > u for u < 1 and G{u) < u for w > 1. In 
the numerical simulations [25] one of the simplest choices for G was implemented 
as 



G{u) = 



(u + l)/2 foru>0 
{u - l)/2 for u < 0, 



(10) 



and random initial and periodic boundary conditions were applied. 

In Fig. 4 we show the time evolution of the model for 77 = 2.0. In a short 
time the system reaches an ordered state, characterized by a spontaneous broken 
symmetry and clustering of the particles. In contrast, for -q = 4.0 the system 
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Fig. 4. The first 3000 time steps of the Id SPP model [L = 300, N = 600, r] = 2.0 (a) 
and rj = 4.0 (b)]. The darker gray scale represents higher particle density. 
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Fig. 5. Phase diagram in the p — rj plane, the critical line follows rjc(p) ^ p'^ . The solid 
curve represents a fit with k = 1/4. 
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would remain in a disordered state. The p — rj phase diagram is shown in Fig. 5. 
The critical line, ric{g), follows (5) with k w 1/4. 

6 Analytical studies of SPP systems 

To understand the phase transitions observed in the models described in the 
previous section, efforts has been made to set up a consistent continuum theory 
in terms of v and p, representing the coarse-grained velocity and density fields, 
respectively. The first approach [2f] has been made by J. Toner and Y. Tu who 
investigated the following set of equations 

dtv + (vV)v = av - /3|vpv - VP + Dl\7{Vv) + DiV^v + D2ivV)^v + ^ 
9tP + V(pv)=0, (11) 

where the a,P > terms make v have a nonzero magnitude, I?l,i,2 are diffu- 
sion constants and ^ is an uncorrelated Gaussian random noise. The pressure P 
depends on the local density only, as given by the expansion 

pip) = j2<j^{p-por- (12) 

n 

The non-intuitive terms in Eq.(ll) were generated by the renormalization pro- 
cess. 

Tu and Toner were able to treat the problem analytically and show the 
existence of an ordered phase in 2d, and also extracted exponents characterizing 
the density-density correlation function. They showed that the upper critical 
dimension for their model is 4, and the theory docs not allow an ordered phase 
in Id. 

However, as we have shown in the previous section, there exist SPP systems 
in one dimension which exhibit an ordered phase for low noise level. Such sys- 
tems can not belong to the universality class investigated in [21]. This finding 
motivated the construction of an other continuum model, which can be derived 
from the master equation of the Id SPP system [25]: 

d,U - f{U) + p^dlU + a^M}^Ml. + Q^ (13) 

p 

dtp ^ -vod^ipU) + Dd^p, (14) 

where U{x, t) is the coarse-grained dimensionless velocity field, the self-propulsion 
term, f{U), is an antisymmetric function with f{U) > for < U < 1 and 
/([/) < for U > 1. The noise, (, has zero mean and a density-dependent stan- 
dard deviation: C,'^ — a^ /pr^. These equations are different both from the equilib- 
rium field theories and from the nonequilibrium system defined through Eqs(ll). 
The main difference comes from the nature of the coupling term {dxU){dxp)/p 
which can be derived from the microscopic alignment rule (1) [26]. Note, that the 
noise also has different statistical properties from the one considered in (11). For 
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a = the dynamics of the velocity field U is independent of p and with an ap- 
propriate choice of/ Eq.(13) becomes equivalent to the <P^ model describing spin 
chains, where domains of opposite magnetization develop at finite temperatures 

[!]• 

7 Linear Stability Analysis 

To study the effect of the nonlinear coupling term a{dxU){dxp)/ p, we now in- 
vestigate the development of the ordered phase in the deterministic case (a = 0) 
when c, D <S^ 1 holds. It can be seen that stationary domain wall solutions exist 
for any a. In particular, let us denote by p* and U* the stationary solutions 
which satisfy p*i±oo) = 0, U*{±oo) = Tl, U*{x) < for a; > and U*{x) > 
for X < 0. These functions are determined as 

ln^ = ^ rU*{x')dx' (15) 



P*(0) D Jo 

and 

^^dlU* ^-,f{U*)-a^U*d,U*. (16) 

Although stationary solutions exist for any a, they are not always stable against 
long wavelength perturbations, as the following linear stability analysis shows. 

Let us assume that at t = we superimpose an u{x, t = Q) perturbation over 
the U*{x) stationary solution. Since c, I? <C 1 the dynamics of p is slow, thus 
p{x,t) = p*{x) is assumed. The stationary solutions are metastable in the sense 
that small perturbations can transform them into equivalent, shifted solutions. 
Thus here by linear stability or instability we mean the existence or inexistence 
of a stationary solution to which the system converges during the response to a 
small perturbation. To handle the set of possible metastable solutions we write 
the perturbation u in the form of u as 

U{x,t) ^U*{x)+u{x,t) ^U*{x-^{t))+u{x,t), (17) 

i.e., 

u{x,t)^u{x,t)+^{t)a{x), (18) 

where a < denotes dxU* and the position of the domain wall, ^{t), is defined 
by the implicit 

u{m,t)^o (19) 

equation. As t/*(0) — 0, from (19) we have u{^,t) = 0. The usage of ^ and u is 
convenient, since the stability of the stationary solution U* is equivalent with 
the convergence of ^(t) as t ^ oo. 

Substituting (17) into (14) and taking into account the stationarity of U* we 
get 

-a{x ~ ^)i + dtuix) ^ (/' o U*){x - ^)u{x) + p^dlu{x) + a^h{x - ^), (20) 
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with h — adi Inp* > 0. To simplify Eq.(20 let us write (/' o U*){x) in the form 
of g{x) - goo-, where g^o = -limx-,±oof'{U*{x)) = -/'(±1). Furthermore, a 
moving frame of reference y = x — ^ and the new variable v{y) = u{x) can be 
defined. With these new notations Eq.(20) becomes 



-a{y)i + dtv{y) = g{y)v{y) - goov{y) + n'^dlv{y) + a£,h{y). 



(21) 



The time development of ^ is determined by the u(^) — v{0) — condition 
yielding 

-a(0)e = tJ-'^dlv{0) + a^h{0). (22) 

By the Fourier-transformation of Eq. (21), treating the gv term as a pertur- 
bation, for the time derivatives of S.{t) and the n-th Fourier moments Vn{t) ~ 
9"w(0, t) one can obtain [27] 



(23) 



Expression (23) can be further simplified using the relations cim <^ «« and 
hm ^ ft-n for TO > n and the approximate solutions for the A growth rate of the 
original u perturbation we found [27] 
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(25) 
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If a = 0, then A+ = and A_ = gca^o — /i^a2 < 0. However, for certain a > 
values A+ > can hold obtained as either g < or 6 > 0, i.e., 



or 



a > ad ^ D 



d > ac.2 — D 



3.goo - 2,g(0) 
2ca(0) 

250c - g(0) 



(28) 



(29) 



a(0)(c-D)' 

respectively. Thus for a > ac — min{ac^i, 010,2) the stability of the domain wall 
solution disappears. 

The instability of the domain wall solutions gives rise to the ordering of the 
system as the following simplified picture shows. A small domain of (left mov- 
ing) particles moving opposite to the surrounding (right moving) ones is hound 
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to interact with more and more right moving particles and, as a result, the do- 
main wall assumes a specific structure which is characterized by a buildup of the 
right moving particles on its left side, while no more than the originally present 
left moving particles remain on the right side of the domain wall. This pro- 
cess "transforms" non-local information (the size of the corresponding domains) 
into a local asymmetry of the particle density which, through the instability of 
the walls, results in a leftward motion of the domain wall, and consequently, 
eliminates the smaller domain. 

This can be demonstrated schematically as 

A B 

where by > (<) we denoted the right (left) moving particles. In contrast to 
the Ising model the A and B walls are very different and have nonequivalent 
properties. In this situation the B wall will break into a Bi and i32, moving in 
opposite directions, Bi moving to the left and B2 moving to the right, leaving 
the area Bi — B2 behind, which is depleted of particles. 

>>>>>>>>>>><<<<< >>>>>>>>>> 

A Bi B2 

At the A boundary the two type of particles slow down, while, due to the insta- 
bility we showed to be present in the system, the wall itself moves in a certain 
direction, most probably to the right. Even in the other, extremely rare case 
(i.e., when the A wall moves to the left), an elimination of the smaller domain 
{A — Bi) takes place since the velocity of the domain wall A is smaller than the 
velocity of the particles in the " bulk" and at Bi where the local average velocity 
is the same as the preferred velocity of the particles. Thus, the particles tend to 
accumulate at the domain wall A, which again, through local interactions leads 
to the elimination of the domain A— Bi. 

It is easy to see that the U — ±1 solutions are absolute stable against small 
perturbations, thus it is plausible to assume that the system converges into those 
solutions even for finite noise. 
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